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O , Abstract 

psj I We study the propagation of two-dimensional finite-amplitude shear 

' waves in a nonlinear pre-strained incompressible solid, and derive sev- 

^ ■ eral asymptotic amplitude equations in a simple, consistent, and rig- 

^ , orous manner. The scalar Zabolotskaya (Z) equation is shown to be 

• i— I ' the asymptotic limit of the equations of motion for all elastic gen- 

^ ■ eralized neo-Hookean solids (with strain energy depending only on 

, the first principal invariant of Cauchy-Green strain). However, we 

show that the Z equation cannot be a scalar equation for the propaga- 
tion of two-dimensional shear waves in general elastic materials (with 
strain energy depending on the first and second principal invariants 
of strain). Then we introduce dispersive and dissipative terms to de- 
duce the scalar Kadomtsev-Petviashvili (KP), Zabolotskaya-Khokhlov 
(ZK) and Khokhlov-Zabolotskaya-Kuznetsov (KZK) equations of in- 
compressible solid mechanics. 
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1 Introduction 



The Zabolotskaya (Z) equation describes the propagation of finite-amphtude, 
two-dimensional, hnearly-polarized, shear waves in nonhnear sohds. It can 
be presented in the following form (due to Cramer and Andrews (2003)): 

{Ur + f3U^U^)^ + lU,, = 0, (1.1) 

where U = f/(x, t) is a measure of the amplitude of the wave, x is a spatial 
variable measured in a frame moving with the wave, ?7 is a spatial variable 
transverse to the wave propagation, r is a temporal variable, and P charac- 
terizes the material properties of the solid. Zabolotskaya (1986) derived it 
via a perturbation scheme from the general equations of elastodynamics. It 
is a counterpart to the Zabolotskaya-Khokhlov (ZK) equation, itself describ- 
ing the propagation of weakly nonlinear, weakly diffracting, two-dimensional 
sound beams (Zabolotskaya and Khokhlov, 1969): 

{Ur + PUU^)^ + lU,, = 0. (1.2) 

(Note that this form is also due to Cramer and Andrews (2003), and that 
the nonlinear parameters /3 are not the same in each equation). The history 
of the ZK equation and its variants has recently been retraced by Rudenko 
(2010). 

There are two major differences between the Z and ZK equations. One 
is obviously the cubic term in fll.ip versus the quadratic term in (11. 2p . The 
other is that the ZK equation governs the motion of the longitudinal displace- 
ment and is a scalar equation, while the Z equation governs the propagation of 
plane-polarized shear waves and is thus a vector equation a priori. Nonethe- 
less, recent papers (Enflo et al., 2006; Wochner et al., 2008; Pinton et al., 
2010) have investigated the propagation of linearly-polarized nonlinear shear 
waves in incompressible solids. 

Here we aim at finding what restrictions must be imposed, if at all, on the 
form of the strain energy W (say), in order to obtain a scalar Z equation. We 
also include the effects of pre-strain, similar to Cramer and Andrews (2003) 
in the compressible case. We point out that another problem of interest re- 
lated to these model equations concerns the mathematical validation of the 
various approximations. This problem has been solved in the framework of 
the Navier-Stokes and Euler system by Rozanova- Pier rat (2009). Moreover, 
because we are mainly interested in the understanding the physical under- 
pinnings at the basement of various modelling assumptions, we adopt here 
a direct method for the derivation of the model equations, not a general re- 
ductive perturbation scheme such as the one proposed by Taniuti (1990) for 
example. 
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In ^we establish that a scalar Z equation can indeed be generated for 
the whole class of the so-called generalized neo-Hookean solids, for which the 
strain-energy density does not depend on the second principal invariant of 
the Cauchy-Green strain tensor. We then try to generalize the study to the 
case of a generic fourth- order incompressible solid, for which W is expanded 
up to the fourth order in invariants of the Green strain tensor. We show 
in §dB.2l that there, the Z equation is necessarily a vector equation that is, 
only plane-polarized — not linearly-polarized — two-dimensional shear waves 
are possible. 

When dissipative and dispersive effects are taken into account for the 
propagation of longitudinal waves, the ZK equation is generalized to 

[Ur + (3UU^ + P'UUrr + f^'U^rr)^ + \Ur,r, = 0, (1.3) 

where /?' and /3" are other material constants. This scalar equation covers 
the purely dissipative (at /3" = 0) Khokhlov-Zabolotskaya-Kuznetsov (KZK) 
equation (Kuznetsov, 1971) and the purely dispersive Kadomtsev-Petviash- 
vili (KP) equation (Kadomtsev and Petviashvili, 1970). In §4] we consider 
dissipative and dispersive effects in generalized neo-Hookean solids. We de- 
rive the counterpart to (11. 3p for shear waves, as 

{Ur + PU^U^ + P'UUrr + P"Urrr)^ + \U^r, = 0, (1.4) 

a scalar equation which covers the Z equation, the modified Khokhlov-Za- 
bolotskaya-Kuznetsov (mKZK) equation and the modified Kadomtsev-Pet- 
viashvili (mKP) equation. 

The aim of this paper is to provide a simple derivation of these scalar 
equations for incompressible materials. This quest is largely motivated by 
the pressing need to model acoustic wave propagation in soft tissues and 
gels, where the combined effects of incompressibility, pre-strain, viscosity, 
and dispersion cannot be ignored. This is a topic of fundamental relevance 
to ultrasound-based techniques of medical imaging such as elastography. 

2 Basic Equations 

Consider a body of incompressible isotropic elastic material, initially at rest 
in a reference configuration Bq, say. For the study of anti-plane shear mo- 
tions, we choose a system of Cartesian coordinates {Xi,X2,X^) to denote 
the position of a particle in the body in Bq. We then call {xi,X2,xz) the 
Cartesian coordinates, aligned with {Xi,X2,Xz), corresponding to the cur- 
rent position X in the deformed configuration B, say. We study the anti-plane 
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motion 

Xi = A-l/2Xi, X2 = A-1/2X2, X3 = AX3+M(Xi,X2,t), (2.1) 

where A is a positive constant — the axial pre-stretch — and n is a three-times 
continuously differentiable function — the shear motion. 

We call F = dx/dX the associated deformation gradient, and B = 
FF^ the left Cauchy-Green strain tensor. The first principal isotropic strain 
invariant, Ji = iiB, is here 

h = \^ + 2\~^ + ul + ul, (2.2) 

and the second principal isotropic strain invariant, I2 = ^[^i — tr(-B^)], is 
equal to 

l2 = + 2\ + \~\ul + ul), (2.3) 

where ui = du/dXi, U2 = du/dX2. In all generality, the strain energy 
density is a function of Ji and I2 only, and the Cauchy stress tensor cr is 
derived as (see e.g. Ogden, 1984), 

(T = -pi + 2WiB -2W2B-^, (2.4) 

where Wi = dW/dli, W2 = dW/dl2, and j9 is a Lagrange multiplier intro- 
duced by the constraint of incompressibility (which imposes that det F = 1 
at all times). 

The equations of motion, in absence of body forces, are 

div a = pd'^x/dt^, (2.5) 

where p is the mass density, which remains constant due to the incompress- 
ibility constraint. Here, the equations of motion reduce to the following 
system of partial differential equations (see Knowles (1976) for the static 
case) , 

qj + 2{W2UjUi). = 0, J = 1,2 (2.6) 
2[{Wi + A-^1^2) Ui]^ = pd\/dt^, (2.7) 

where summation over the repeated index i = 1, 2 is assumed. The auxiliary 
quantity q is given by 

q = Xp{Xi, X2) - 2Wi - 2(A2 + A"^ + ul + ul)W2. (2.8) 
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3 The Z equation 



3.1 Results for generalized neo-Hookean materials 

The system (12.61) consists of three partial differential equations for the two 
unknown functions p and u. Therefore the system is overdetermined and it 
is not possible, in general, to find a general solution for arbitrary choices of 
the strain energy W. Progress is possible for some special forms of the strain 
energy density, in particular for the Mooney-Rivlin form, which is linear in Ii 
and I2 (Adkins, 1954), but only for static solutions, see Knowles (1976) for 
details. The dynamic situation for antiplane motions is much more complex, 
as discussed by Hayes and Saccomandi (2004) in the Mooney-Rivlin case. 

However, if we consider materials for which W = W{Ii — 3) only — the 
so-called generalized neo-Hookean materials — , then the system (12. 6p reduces 
to a well- determined system 

gi = g2 = 0, 2{W'u,), = pd\/dt^, (3.1) 

where W is the derivative of W with respect to its argument (Ji — 3). 

3.1.1 Derivation of the Z equation 

To derive the Z equation, we assume that W is at least of class C^. We 
assume that the primary propagation direction is in the Xi-direction and 
introduce the following new displacement and new variables 

V = e~^u, X = e^-^i! r = t ± aXi, rj = e^faX^-, (3.2) 

where e is a small, non-dimensional parameter, giving a measure of the mo- 
tion amplitude, the plus or minus sign corresponds to left or right running 
waves, and a is a parameter to be adjusted to cancel all terms to order e (Note 
that the scalings used here mirror those of Zabolostkaya (1986), and give a 
consistent asymptotic expansion). To obtain the Z equation^ we expand the 
equation (13. ip . that is 

2(iy'un + W[ux + W'ux2 + W^U2) = pd'^u/dt^, (3.3) 

to order (note again that the subscript ( )j here denotes derivation w.r.t. 
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First, we compute the partial derivatives as follows, 
d . d 2 ^ 



±a— — h e 



dXi dr dx 

d 



8X2 ^ ^87]' 
W ~ 8^' 

8^ (9^ 

(3-4) 



8X1 " "8r] 
which leads to 

Ml = ±eaVr + e^i'x' '^2 = e^y/aVrj, 

uii = ec?VrT ± Se'^f^j.r + O (e^) , ^22 = e'^av^jr^. (3.5) 
Second, we expand W in e. Since 

/, - 3 = A2 + - 3 + e2a\;2 + o (^4) ^ (g g) 

we have 

H^' = u;(i)+eV2)ay + 0(e^) (3.7) 

where w^^) = W'{X^ + 2A-i - 3) and w^^) = vr"(A2 + 2\~^ - 3). Third, we 
use the expansion for W and Ui to compute each term in (13. 3p up to order 
e^. That is 

Wi = ±eV^a' {vl)^ + O {e') , = e'w^'^a'/' {v'^)^ + O (e^) . (3.8) 

Equation (13. 3p reads now 

evrr{'^w'^^^a^ - p) + 2e^a [w^^\v^^ ± 2v^r) + w^'^'^ a^{vl)r] + O (e^) = 0.(3.9) 



Terms to order e disappear by choosing = p/{2w^^^), which gives to order 

Vr,^±2v^r + ^(y''{vl)r = 0. (3.10) 



w 



Finally, taking a derivative of this equation with respect to r and using 
the substitution U = v^. leads to the scalar Z equation: 

{±U^ + PU^Ur)^ + lU^rj = 0, (3.11) 
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where /3 = 3u;''^^a^/p. 

For example consider the Yeoh strain energy described by 

W = (/i/2) [ih - 3) + 7(/i - 3)V2] , (3.12) 

where /i > is the initial shear modulus and 7 > is a nonlinear elastic 
constant. In this case we find 

= Ml + 7(A^ + 2A--3)]' ^ = V"'- 

Note that the apparent difference between (13. lip and the Z equation (II. ip 
is simply due to the choice of independent variables. Here, we chose the same 
independent variables as Norris (1998), where r is a retarded/advanced time, 
X is the spatial variable in the propagation direction, and 77 is the transverse 
spatial variable. 



3.1.2 Special solutions 

Following Gates and Crighton (1990) and Sionoid and Gates (1994), we per- 
form the transformation 

f/(x,r,r/)=G(x,C), where ( = ^ + , (3.14) 

which reduces the Z equation (13.111) to an equation which, once integrated 
in (, yields 

G^ + G^Gc + —G = 0. (3.15) 

The solution of this one- dimensional equation does not provide a general so- 
lution to the two-dimensional equation (13. lip ; rather, it gives the solution to 
the particular Cauchy problems where initial data are given on the parabolic 
surfaces C = const 

A further change of function (from G to G) and of variables (from x to 

x), 

G{xX) = x'^'G{xX), X = lnx, (3.16) 
reduces equation (13. lip to the autonomous inviscid generalized Burgers equa- 
tion 

G^ + G^G^ = 0. (3.17) 

The general integral of this equation is G = F (^(^ — G^x^, where F is an 

arbitrary function. This allows us to recover the following class of exact 
implicit solution for the Z equation (13. lip . 



U = x-"^F (r+^- /3[/^xe^) • (3.18) 



It is quite remarkable that the transformation fl3.14p should work for the 
Z equation (cubic nonlinearity) , because it was initially intended for the KZK 
equation (quadratic nonlinearity). Through it, a huge class of exact solutions 
is generated for the Z equation. 

General methods to obtain exact solutions for such equations include 
reductions methods such as the one based on Lie group analysis and its 
generalizations. Those methods have been systematically applied to the ZK 
equations with quadratic nonlinearities (see for example Bruzon et al. (2009) 
and Tajiri (1995)), but it seems that there are very few results regarding the 
case of a cubic nonlinearity, of interest to incompressible elastic materials. 



3.2 Results for general incompressible solids 



We return to an important point evoked in the Introduction: that the prop- 
agation of finite-amplitude transverse waves cannot be governed by a scalar 
Z equation (i.e., linearly-polarized shear waves exist only in special solids, 
such as the generalized neo-Hookean materials). 

To show this, we take solids with the following Rivlin strain energy. 



W = Cwih - 3) + Coiih - 3) + C2o(/i - 3)2, 



(3.19) 



(where Cio, Cqi, C20 are constants), which covers the most general model of 
fourth-order incompressible elasticity (Ogden, 1974), 



W = fitT [E^] + 4 tr {E^) + D {tT{E')Y 



(3.20) 



with the identifications (Destrade et al., 2010): 
/i = 2(Cio + Coi), A = -8(Cio + 2Coi), D 



2(G 
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3Coi + 2C2o)- 
(3.21) 

Here E = {F^F — I)/2 is the Green strain tensor and fi, A, and D are the 
second-, third-, and fourth-order elasticity constants, respectively. 

With this strain energy, the first and second equations of (12.61) read 



+ 2Coi 
qx, + 2Coi 



u 



u 



X2 



X2 



+ ( 



0, 



0, 



(3.22) 



where q = q/L and Xi, X2 are defined as in the previous sections. 

Guided by the results of the previous subsection, we know that we must 
perform the changes of function and variables defined by (13. 2|) . if we aim at 
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deriving the Z equation from the last equation in f l2.6p . Here, those changes 
give the first terms in the expansion of the two equations fl3.22p as 



e^q^ ± aqr ± 2Coie^a^ (U^)^ + . . . = 0, eg^ + 2Coie^a {U^Ur^^ + . . . = 0. 

(3.23) 

Now the first of these tells us that q must be of order e^: q = e^Q, say. Then 
the equations f l3.22p reduce further to 

Qr + 2Coia (U^) ^ = 0, Q, + 2Coia (ifrU,) ^ = 0, (3.24) 

at order e^. Clearly these expressions for the derivatives of Q are not com- 
patible. 

We may thus conclude that linearly-polarized waves are not possible in 
general incompressible elastic materials. Indeed, Horgan and Saccomandi 
(2003) explain that when an anti-plane shear deformation occurs in a general 
incompressible material, it is necessarily accompanied by secondary motions 
associated with in-plane deformations. This coupling increases dramatically 
the complexity of the calculations and is an obstruction to the propagation 
of linearly polarized two-dimensional shear waves. In general incompressible 
materials it is necessary to take into account the in-plane motions coupled 
to the two-dimensional shear wave. 

In fact, Wochner et al. (2008) acknowledge as much when they correct 
Enflo et al. (2006) and their study of linearly-polarized nonlinear shear 
waves, and conclude that "Simply stated, quadratic nonlinearity prevents 
the existence of a pure linearly polarized shear wave beam." However, they 
then state that "Cubic nonlinearity can be more pronounced than quadratic 
nonlinearity due to the quasiplanar nature of wavefronts in directional sound 
beams", and go on to "consider only the effects of cubic nonlinearity." But 
ignoring quadratic terms in favour of cubic terms runs contrary to the very 
idea of asymptotic expansions. Upon examination of their equation of motion 
(11), we see that the only way to make quadratic terms vanish is to take 
4fi + A = 0. According to (13.210 here, this is equivalent to saying that 
Cqi = 0, making W in fl3.20p belong to the class of generalized neo-Hookean 
materials. 

4 Dissipation and Dispersion 

In our treatment of dissipative and dispersive effects, we use the following 
constitutive law for the Cauchy stress: 

o- = a^ + a^, (4.1) 
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where 



cr^ = + /i [1 + 7(/i - 3)] ct"" = uAi + aiiA2- Al), (4.2) 

which separates the purely elastic part cr^ of the stress from a viscoelastic 
part (T^. 

More specifically, cr^ coincides with the Cauchy stress of an elastic Yeoh 
material with strain energy (I3.12p . We make this choice for convenience, but 
it is clear that once the asymptotic expansions are performed below, the Yeoh 
form recovers the entire class of generalized neo-Hookean solids, through the 
identifications: 

+ 7(A2 + 2A-^ - 3)] = 2w;(^), ji'y = 2w^'^\ (4.3) 

The viscoelastic part of the stress, cr^, is the stress of a second-grade 
viscoelastic solid, where u > 0, ai > are constants and Ai, A2 are the first 
two Rivlin-Ericksen tensors, defined by 

Ai = L + L^, A2 = Ai + AiL + L^A,, (4.4) 

with L = FF^^, the displacement gradient. Here the superposed dot rep- 
resents the material time derivative. The first part of cr^ is uAi = 2uD, 
where D is the stretch tensor; that part describes dissipation by Newtonian 
viscosity effects, and coincides with the dissipative term in the Navier-Stokes 
theory. The second part of cr^ is ai{A2 — A^); it describes dispersion by 
accounting for the intrinsic characteristic lengths of the solid, see Destrade 
and Saccomandi (2006) for details. 

Note that the temporal evolution of the canonical energy E{t) for our 
sohd is given by E{t) = -2 /^^ iy{D ■ D)dV < 0, where fi* is 

a given region occupied by the material at time t. It is independent of ai 
(Fosdick and Yu, 1996). 

With respect to two-dimensional anti-plane motion, the derivation of the 
governing equation is not overly complicated by the additional terms. For 
ease of exposition, we consider that the solid is not pre-strained: A = 1 (It 
follows that once the asymptotic expansions are performed, the identifications 
(Oj) reduce to /i = 2w^^\ 7 = 2w^^'> /w^^\) 

We find that the equations = (a = 1,2) still apply here, whilst the 
last equation in (13. ip is changed to 

[fiup + fij{ul + ul)ui3 + uupt + aiuptt] p = putt- (4.5) 

In dimensionless form it becomes 

Uii + 7[(mi + U2)ui\i + -^Uni + ^WiitT - ^"it = 0' (4.6) 
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using the usual scalings. 

When the sohd is not dissipative (z/ = 0) and not dispersive (ai = 0), the 
changes of function and variables (13. 2p give a governing equation at order e^, 
when a is chosen as a = ^pL^inT'^). Here the expansions of the dissipative 
and dispersive terms in (14. 6 p start as 



Q^i - N Oil f pL"^ . ,^ 

(^11« + «22ti) = -7^ f/rrrr + • • • • (4.7) 



Clearly, in order to be incorporated into the governing equation for the mo- 
tion, the dissipative and the dispersive effects must be at least of order e^: 
V = e^v + . . . and ai = e^di + . . . say, where z> and di are constants. With 
that assumption, we find the following form of the equation of motion for the 
velocity U = Ur, 



where 



[±U^ + PU'^Ur + P'Urr + P"Urrr)^ + \Ur,, = 0, (4.8) 



3 1 1 



(4.9) 

This equation covers not only the Z equation (at /?' = (5" = 0), but also the 
dissipative mZK equation (at /3" = 0) and the dispersive mKP equation (at 
/3' = 0). 

To the best of our knowledge this the first time that a mKP equation is 
derived in the continuum mechanics of homogeneous solids. 
If we consider the diffusive version of (14. 8p . i.e. 

{±U^ + f3U^Ur + f3'Urr)^ + \Ur,^ = 0, (4.10) 

and we seek solutions which propagate at an angle 7 with respect to the r 
direction, it is possible to obtain an interesting reduction. To this end let us 
consider solutions which depend only on x and f = r+ \/2ri tan(7) — x tan2(7) 
to reduce (I4.10p to 

{±U^ + f3U^Uf + f3'Uff)^ = 0. (4.11) 

This is the Nariboli-Lin modified diffusive Burgers' equation (Nariboli and 
Lin, 1973). The plane wave solution of this equation is computed as 



U-' = ^ + exp 



[T - Vox) 



11 



where Vq is the speed of the plane wave. This solution is quite different from 
the plane wave solution of the classical Burgers equation which is expressed 
in terms of the hyperbolic tangent function. 

If we consider the purely dispersive case and we rescale some of the vari- 
ables and U = XiU, x = and 77 = Xsfj it is always possible to choose the 
scaling such that 



Ufjfj 



0, 



(4.12) 



the modified Kadmotsev-Petviashvili equation. A connection between this 
equation and the Modified Korteweg de Vries equation may be established 
with a similar transformation to the one used to reduce (14.101) to (14. lip . 
It is also possible to seek direct solutions of (14.121) in the form 

IJ = F{Q, C = af + bfj ~ ux- 

This ansatz produces the reduction 



a'F" 



(4.13) 



where Cq, Ci are integration constants. 

If we consider solitary waves solutions the asymptotic boundary condi- 
tions require Cq = Ci = and the usual energy integral is given by 

2 T2 



1 IP 



C2, 



where again C2 is an integration constant. It is clear that it must be C2 
and if acj — 6^ > an explicit solitary wave solution is given by 

1/2 







F(C) = {Q{au- b'^) sech^ 



au 



62 



2a 



(C - Co) 



where (q is an integration constant. 

If we consider in (I4.13P cq = and ci 7^ it is also possible to find cnoidal 
wave solutions in terms of elliptic functions. 



5 Conclusions 

We provide a clear and complete derivation of various model equations for 
shear waves in the framework of nonlinear incompressible elasticity. The ra- 
tionale for investigating this type of problem is dictated by the use of imag- 
ing techniques of improved quality, such as latest generation elastography, to 
evaluate the nonlinear material properties of soft tissues. 
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We show that it is possible to derive a pure (scalar) Z equation only in the 
special case of generalized neo-Hookean materials. For general incompressible 
materials the Z equation is coupled with a non-negligeable transverse motion, 
an important issue which is often overlooked in the literature. Clearly, a 
rigorous evaluation of the validity of the perturbative expansion requires an 
extension of the existence, uniqueness and stability results that have been 
provided so far only for model equations with quadratic nonlinearity. 

We point out that general structure of the various equations here obtained 
is as follows 

-^^(f/) + cV^f/ = 0, 
dx 

where c is a suitable material parameter, V_l is the second order derivative 
(Laplace operator) in an orthogonal direction to the beam axis and S{U) = 
is the equation governing plane wave propagation (inviscid Burgers, Burgers 
or KdV equation). It would be an interesting problem to study the various 
transformation and group analysis properties of such general equation to find 
exact and explicit solutions for various model equations following a general 
and complete scheme. 
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